# ------------------------------------- #
# Replication code for:
#
# Pomeroy, Caleb "The Damocles Delusion: The Sense of Power Inflates Threat  
# Perception in World Politics," International Organization.
#
# This script reproduces the survey experimental results from China.
# ------------------------------------- #


# --- set your working directory --- #
setwd("~/Downloads/damocles_delusion_replication/")

# --- libraries --- #
library(ggplot2)
library(plyr)
library(texreg)

# --- load data --- #
china_survey <- readRDS("data/china_nuclear_proliferation_data.rds") 
nrow(china_survey) # the China survey consists of N = 880 Beijing-based respondents gathered by Qualtrics
china_survey$nuke_threat <- as.numeric(china_survey$nuke_threat) # threat perception on a ten-point scale (higher = more threatening)
china_survey$nuke_attack <- as.numeric(china_survey$nuke_attack) # support for preventive attacks on a ten-point scale (higher = more support)
china_survey$power_condition <- factor(china_survey$power_condition, levels = c("nuke_control", "nuke_high")) # subjects' symmetric ("control") versus high power ("high") condition assignment
china_survey$aggressive_intentions <- scales::rescale(china_survey$aggressive_intentions, c(0,1)) # seven-point pre-treatment (dis)agreement with the prior belief that other countries harbor aggressive intentions towards China (higher = more aggressive intentions) 
china_survey$nat_attach <- scale(china_survey$nat_attach)[,1] # seven-point disagreement with the notion that "Being Chinese is important to how I feel about myself" (higher = more nationally attached)
china_survey$education <- factor(china_survey$education, levels = c("less_bachelors", "bachelors")) # whether the respondent has completed at least a four-year college degree ("bachelors")  or not ("less_bachelors")
china_survey$gender <- factor(china_survey$gender, levels = c("female", "male")) # self-reported male versus female respondent
# china_survey$age # age is on an eight point scale, where 1 = 18-24, 2 = 25-34, 3 = 35-44, 4 = 45-54, 5 = 55-64, 6 = 65-74, 7 = 75-84, 8 = 85 or older
china_survey$mi_war <- scale(china_survey$mi_war)[,1] # seven-point (dis)agreement with the notion that "War is unfortunate but sometimes the only solution to international problems" (higher = more militant)

# --- sample demographics --- #
# below are sample demographic stats mentioned in the appendix
table(china_survey$gender)/nrow(china_survey) # 52% male
sum(table(china_survey$age)[1:2])/nrow(china_survey) # 35% aged 18-34
sum(table(china_survey$age)[3:5])/nrow(china_survey) # 54% aged 35-64
sum(table(china_survey$age)[6:8])/nrow(china_survey) # 11% aged 65 or older
table(china_survey$education)/nrow(china_survey) # 64% have a four-year degree (or higher)

# --- regression models --- #
# main effects of power on threat perception and support for preventive force
summary(threat_model <- lm(nuke_threat ~ power_condition + gender + age + education + 
                             nat_attach + mi_war, china_survey))
summary(attack_model <- lm(nuke_attack ~ power_condition + gender + age + education + 
                             nat_attach + mi_war, china_survey))

# --- appendix table A1
screenreg(list(threat_model, attack_model))
# texreg(list(threat_model, attack_model),
#        booktabs = T,
#        caption.above = T,
#        caption = "OLS Estimates: China Experiment Results",
#        custom.coef.names = c("(Intercept)",
#                              "High Power Treatment",
#                              "Male",
#                              "Age",
#                              "4-Year Degree",
#                              "National Attachment",
#                              "Militant Orientation"),
#        stars = c(.001, .01, .05, .10),
#        symbol = "\\wedge",
#        custom.model.names = c("Threat Perception",
#                               "Preventive Attacks"))

# --- format for plotting purposes
mod_list <- list("threat_model" = threat_model, 
                 "attack_model" = attack_model)
results_df <- data.frame()
for(i in 1:length(mod_list)){
  
  mod_i <- mod_list[[i]]
  
  results_df <- rbind(results_df,
                      data.frame(ci_95_high = confint(mod_i, level = .95)[,2],
                                 ci_95_low = confint(mod_i, level = .95)[,1],
                                 ci_90_high = confint(mod_i, level = .90)[,2],
                                 ci_90_low = confint(mod_i, level = .90)[,1],
                                 est = coef(mod_i),
                                 variable = names(coef(mod_i)),
                                 model = names(mod_list)[i]))
}
results_df$variable <- revalue(results_df$variable, c("power_conditionnuke_high" = "High Power\nTreatment",
                                                      "gendermale" = "Male",
                                                      "age" = "Age",
                                                      "mi_war" = "Militant\nOrientation",
                                                      "educationbachelors" = "4-Year Degree",
                                                      "nat_attach" = "National\nAttachment")) 
results_df$model <- factor(results_df$model, levels = c("threat_model", "attack_model"))
levels(results_df$model) <- c("Threat Perception", "Preventive Attacks")
results_df <- results_df[!results_df$variable%in%c("(Intercept)", "Age", "Male", "4-Year Degree"),]
results_df$variable <- 
  factor(results_df$variable, rev(c("High Power\nTreatment", "Militant\nOrientation",
                                    "National\nAttachment")))
plot_col <- ifelse(results_df$variable == "High Power\nTreatment", "#54545d", "#99998e")

# --- figure 1 from the main text
ggplot(results_df, aes(x = est, y = variable, color = plot_col)) +
  geom_vline(xintercept = 0, linetype = "solid", color = "black", size = .8) +
  geom_linerange(aes(xmin = ci_95_low, xmax = ci_95_high), size = 1, color = "black") +
  geom_linerange(aes(xmin = ci_90_low, xmax = ci_90_high), size = 2.8) +
  geom_point(size = 5, color = "black") +
  geom_point(size = 4) +
  facet_wrap(~model, nrow = 1) +
  theme_bw() +
  labs(x = "Estimate", y = NULL) +
  scale_color_manual(values = plot_col) + 
  theme(axis.title.x = element_text(size = 18, margin = margin(t=15)),
        axis.text.y = element_text(size = 14, color="black"),
        axis.text.x = element_text(size = 13),
        strip.text = element_text(size = 16, margin = margin(r=8,l=8,t=8,b=8)),
        legend.position = "none",
        panel.spacing = unit(1.2, "lines"),
        panel.background = element_rect(fill = "#f6f6f5"),
        strip.background = element_rect(fill = "gray90", color = "black", size = 1.5),
        plot.margin = margin(r=75, l=55, t= 15))


# --- moderation analysis --- #
# moderation analysis of power treatment x prior beliefs about aggressive intentions
summary(threat_moderation <- lm(nuke_threat ~ power_condition*aggressive_intentions + gender + age + education + 
                             nat_attach + mi_war, china_survey))
summary(attack_moderation <- lm(nuke_attack ~ power_condition*aggressive_intentions + gender + age + education + 
                             nat_attach + mi_war, china_survey))

# --- appendix table A2
screenreg(list(threat_moderation, attack_moderation))
# texreg(list(threat_moderation, attack_moderation),
#        booktabs = T,
#        caption.above = T,
#        caption = "OLS Estimates: China Experiment Moderation Effects",
#        custom.coef.names = c("(Intercept)",
#                              "High Power Treatment",
#                              "Prior Beliefs (Aggressive Intentions)",
#                              "Male",
#                              "Age",
#                              "4-Year Degree",
#                              "National Attachment",
#                              "Militant Orientation",
#                              "High Power Treatment x Prior Beliefs"),
#        stars = c(.001, .01, .05, .10),
#        symbol = "\\wedge",
#        custom.model.names = c("Threat Perception",
#                               "Preventive Attacks"))

